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Abstract. - We derive a set of general dynamical mean- field equations for strongly interacting 
fermionic atoms in arbitrary slowly-varying potential traps. Our derivation generalizes the 
ansatz of the crossover wavefunction to the inhomogeneous case. The equations reduce to a 
time-dependent Gross-Pitaevskii equation on the BEC side of the resonance. We discuss an 
iteration method to solve these mean-field equations, and present the solution for a harmonic 
trap as an illustrating example, which self-consistently verifies the approximations made in our 
derivation. 



The recent experimental advance in ultracold atoms has allowed controlled studies of 
strongly interacting Fermi gases in various types of potential traps, with the interaction 
strength tunable by an external magnetic field through the Feshbach resonance [1]. For homo- 
geneous strongly interacting Fermi gases at zero temperature, the physics is captured by the 
variational crossover wavefunction [2-4], which interpolates the BCS and the BEC theories. 
For strongly interacting atoms in a potential trap, there are currently two main methods to 
deal with the resultant inhomogeneity: one is the local density approximation (LDA) [4,5], 
which neglects the kinetic terms associated with the spatial variation of the order param- 
eters; the other is based on the numerical simulation of the Bogoliubov-De-Gennes (BDG) 
equations [6,7]. Both of these methods have found wide applications recently [4-7], but 
each of them also has its own limitation: the LDA becomes inadequate in cases where varia- 
tions of the order parameters have significant impacts; the BDG equations take into account 
exactly the spatial variation of the order-parameter, but its numerical solution is typically 
time-consuming, which limits its applications only to very special types of potentials. 

In this work, we develop a different method to describe both the static properties and 
the dynamics of strongly interacting fermionic atoms in arbitrary but slowly varying poten- 
tial traps. Our starting point is a variational state which is a natural generalization of the 
crossover wavefunction to the inhomogeneous and dynamical cases. The key simplification in 
our derivation comes from the assumption that the spatial variation of the order-parameter is 
small within the size of the Cooper pairs. This assumption is similar to the one in the deriva- 
tion of the Ginzberg-Landau equation for the weakly interacting fermions [8] , but we avoid the 
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use of perturbation expansions so that the order parameter here in general does not need to be 
small [8,9]. With such an assumption, we derive a set of dynamical mean-field equations for 
the bare molecular condensate and the Cooper-pair wavefunctions. This set of equations can 
be solved iteratively, and its zeroth-order approximation, which neglects the order-parameter 
variation, gives the LDA result. On the BEC side of the Fcshbach resonance (with the chem- 
ical potential fj, < 0), these mean-field equations can be reduced to a generalized dynamical 
Gross-Pitaevskii (GP) equation [8-10], with the effective nonlinear interaction for the bare 
molecules derived from a fermionic gap equation. When one goes deeper into the BEC region, 
the nonlinear interaction resumes the conventional GP form, and one can derive an effective 
scattering length for the bare molecules. We solve the dynamical mean-field equations for a 
harmonic trap as a simple illustrating example to self-consistently verify the approximations 
made in our derivation. Another recent work also addresses the dynamics of a trapped Fermi 
gas across a Feshbach resonance [11]. The set of equations derived therein are semiclassical 
hydrodynamic equations, which, after linearization, can be applied to calculate the dynam- 
ical properties of the system. In contrast, our work follows the time-dependent variational 
approach. By reducing the dynamical mean field equations to a time-dependent non-linear 
GP equation on the BEC side of the resonance, we provide a complementary perspective to 
the problem. 

Our starting point is the two-channel field Hamiltonian [3, 4, 12] 

H = 12 y^[-V7(2m) + V>,t)]^d 3 r + J ^[-V 2 / (4m) + 7 + 2V{r, t)]V b d 3 r 

+ a J[^ l ^ ] d 3 r + h.c.}+U j *{*{* x * T d 3 r, (1) 

which describes the interaction between the fermionic atom fields If J. (<r =f, j labels atomic 
internal states) in the open channel and the bosonic bare molecule field ^\ in the closed 
channel. In this Hamiltonian, m is the atom mass, and V(r, t) is the trap potential which 
could vary both in space and in time. Note that we have assumed that the trap frequencies for 
a composite boson and for a single atom are the same, so that the potential that a boson feels 
is twice as a single atom does. The bare atom-molecule coupling rate a, the bare background 
scattering rate U, and the bare energy detuning of the closed channel molecular level relative 
to the threshold of the two- atom continuum 7 arc connected with the physical ones a p ,U p , 7 P 
through the standard renormalization relations [4]. The values of the physical parameters 
a P>U p ,"fp are determined respectively from the resonance width, the background scattering 
length, and the magnetic field detuning relative to the Feshbach resonance point (see, e.g., the 
explicit expressions in Ref. [13]). Note that following the standard two-channel model, direct 
collisions between the bosonic bare molecules arc neglected, as their contribution is negligible 
near a broad Feshbach resonance [3,4, 12]. 

At almost zero temperature and with a slowly varying potential V(r,t), the state of the 
Hamiltonian (1) can be assumed to evolve according to the following variational form: 



|$(t)) = Afexp 



J f(r,r',t)^\(r)^l(r')d 3 rd 3 r' exp J b (r, t)*\ (r)d 3 r 



\vac), (2) 



where J\f is the normalization factor, (j)i,(r,t) is the condensate wavefunction for the bare 
molecules, and /(r,r',i) is the Cooper-pair wavefunction. This variational state is a natural 
generalization of the crossover wavefunction to the inhomogeneous and dynamical cases [14]. 
Without the fermionic field, this variational state would have the same form as the one in the 
derivation of the dynamical GP equation for the weakly interacting bosons [10]. 
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To derive the evolution equations for the wavefunctions (f>b(r,t) and f(r,r',t), we follow 
the standard variational procedure to minimize the action S — f dt[(<b\&) — ($|<f>)]/(2i) — 
where $} and H are specified in Eqs. (1) and (2). Under the ansatz state (2), the 
Wick's theorem implies the decomposition [/(*|*{*x* T ) « C/(*|*J)(*j.* T ) (the additional 
Hartree-Fock terms, which only slightly modify the effective V(r, t), are not important when 
there is pairing instability [8] and are thus neglected here). It turns out that to get the 
expression of the action S, the critical part is the calculation of the pair function F* (ri , r 2 , t) = 
(<I'j(r 1 ) , I'j(r 2 )). Under the ansatz state (2), the pair function satisfies the following integral 
equation (we drop the time variables in F* (r x , r 2 , t) and /* (r, r', t) when there is no confusion) 

^*(ri,r 2 ) = .r(r 1 ,r 2 )-y.r(r 1 ,r3).r(r 4 ,r 2 )^(r4,r 3 )(i 3 r3d 3 r4. (3) 

To solve this integral equation, we write both F*(ri,r 2 ) and /*(i"i,r 2 ) in terms of the new 
coordinates r = (ri + r 2 )/2 and r = ri r 2 . Then, we take the Fourier transformation of 
Eq. (3) and its conjugate with respect to the relative coordinate r_. The Fourier transforms 
of F(r, r_) and /(r, r_) are denoted by Fk(r) and /k(r), respectively. In this Fourier trans- 
formation, we assume |<9 r /| <C |3 r _/| and \d r F\ -C |d r _.F|. Physically, it corresponds to the 
assumption that the order parameter is slowly varying within the size of the Cooper pairs. 
Under this assumption, we derive from Eq. (3) and its conjugate the following simple relation 
between the Fourier components 

F k (r) = / k (r)/[l + |/ k (r)| 2 ]. (4) 

This relation is critical for the explicit calculation of the action S. 

We can now express the action S in terms of the variational wavefunctions /k(r,i) and 
for 06 (r, t). From the functional derivatives 5S/8f£(r, t) = and SS/S(f>l(r, t) — 0, we get the 
following evolution equations for /k(r, i) and for <f>b(r, t): 

idtfk = [2e k + flo(r,t)]/k + A(r,t)-A*(r,t)/2, (5) 
id t cf> b = [j + H (r,t)]cf, b + (a/U)A f (r,t), (6) 

where H (r,t) = -V 2 / (4m) + 2V(r,i), £k = h 2 k 2 /2m (m is the atomic mass), A/(r,t) = 
(U/8tt 3 ) J d 3 k/ k (r,t)/ [l + |./k(r,t)| 2 ], and A(r,i) = a<j> b (r,t) + A/(r,i). The two equations 
(5) and (6) represent a central result of this work: they completely determine the evolution of 
the wavefunctions /k and for cj> b , just as the GP equation determines the condensate evolution 
for the weakly interacting bosons. In the stationary case with a time-independent trap, one 
just needs to replace idtfk and id t (j) b respectively with 2/^/ k and 2^(f) bl where \i is the atom 
chemical potential. 

The evolution equations (5) and (6) are a set of coupled nonlinear differential equations. 
They can be solved through direct numerical simulations (for instance, through the split- 
step method), but as the potential V(r, t) is typically slowly varying both in r and in t, 
the following iterative method may prove to be more efficient. In this case, we expect 
both 4> b (r, t)e l2flt and / k (r,i)e l2M * to be slowly varying in r and t. We can then introduce 
the following effective potentials V e ff(r,t) = 0^ 1 (r,t) [—id t — V 2 / (4m) + 2/z] <j> b (r,t)/2 and 
U^ 7 (r, t) = ,/ k x (r, t) [-id t - V 2 / (4m) + 2/z] / k (r, t)/2, both of which should be small. With 
these introduced potentials, we can solve /k(r, t) from Eq. (5) as 



/ k = -(^k-(e k -^ // ))/A 



(7) 
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where Me// = fj, - V(r,t) - V* ff (r,t), Me// = p, - V(r,t) - V eff (r,t), A = a<j) b \l - U(j - 
2fi e ff)/a 2 ], and _E k = (e k - Me//) 2 + ^ 2 - Substituting Eq. (7) into Eq. (6), we get the 
following effective gap equation 

where l/U T = 1/ [U p - a 2 p / ( 7p - 2p eff )] = 1/ [U - a 2 / ( 7 - 2fi eff )] + J d 3 k [1/ (l6. 3 e k )] 
(the latter equality comes from the renormalization relation between j,a,U and 7p , a p , U p [4]). 

Under the zeroth-order approximation, we assume V e ff(r,t) ~ K//( r '*) - °. which leads 
to Me// = Me// = M — ^( r >*) m ^k- In this case, the gap equation (8), together with 
the number equation N = J n(r, t)d 3 r, where N denotes the total atom number and n(r, t) = 
2\4>b\ 2j r (2/87T 3 ) / d 3 k|/k| 2 /(l + |/k| 2 ) is the local atom density, completely solves the problem, 
the result of which corresponds to a solution under the local density approximation in the 
adiabatic limit. Thus we recover the LDA result under the zeroth-order approximation which 
completely neglects V e ff (r, t) and (r, t). It is then evident as how to go beyond the LDA. 

We can use the LDA result ^[°^(r, t), f^\r,t) as the zeroth-order wavefunctions to calculate 
the first-order effective potentials Ujy^r, t) and V^\v, t) through their definition equations. 
Substituting these effective potentials into the gap equation (8) and (7), we can find out the 
next order wavefunctions <f>^ (r, t) , (r, t) . This iterative process should converge if the 
effective potentials are small (i.e., the order parameters are slowly varying in r and t). 

In the following, we consider a different simplification of the basic equations (5) and (6) on 
the BEC side of the resonance with the chemical potential \x < (note that it is not required to 
be in the deep BEC region). On this side, we expect the wavefunctions </>&(r, t) and /k(r, t) to 
have similar dependence on r and t, so we assume V e ff(r, t) ~ V^j(r,t). This approximation 
will be self-consistently tested and we will see that it is well satisfied when fi < 0. Under this 
approximation, /z£ fj = p e f / , while Me/ / can be solved from the gap equation (8) as a function 
of |0f,| 2 . Substituting this solution p e ff into the definition equation of V e ff, we get 

id t <pb = [- V 2 / (4m) + 2V(r, t) + 2p ef /] <f> b . (9) 

This equation has the same form as the dynamical GP equation for the weakly interacting 
bosons except that the collision term is now replaced by a general nonlinear potential Me/ / [10] , 
a function of \<pb\ 2 with its shape determined by the gap equation (8). We have numerically 
solved Eq.(8) for the function /x e //(|</>6| 2 ) at several different detunings for 6 Li, and the shapes 
of these functions are shown in Fig. 1. We can see that Me// becomes almost linear in |0b| 2 
when one goes further into the BEC region. In that limit, Eq. (9) reduces to an exact 
GP equation, and we can define an effective scattering length a e // for the bare molecular 
condensate with cfyie///d(|0b| 2 ) = ^ 7Ta eff/ (2m). This effective scattering length a e // is shown 
in Fig. 1(d) as a function of the field detuning for 6 Li. We should note, however, that 
the effective scattering length for the bare molecules is in general very different from the 
one for the dressed molecules [4,9,15]. The dressed molecules are dominantly composed 
of Cooper pairs of atoms in the open channel (for instance, when the chemical potential 
M ~ 0, the population fraction of the bare molecules is only about 0.1% for 6 Li). As the bare 
molecules have a very low density near the resonance, they in general have a large effective 
scattering length to compensate for that, as is shown in Fig. 1(d). The effective scattering 
lengths for the bare and the dressed molecules coincide with each other only in the deep 
BEC region with the population dominantly in the closed channel. In this limit, we have 



(8) 
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-130.4460 



-130.4462 




Fig. 1 - (a)(b)(c) The effective potential fi e ff as a function of |a</>t,| 2 , with the field detuning B — Bq 
given by — 268G (a), -T07G (b), and 0G (c), respectively. Both u e // and \a</>b\ are in the unit of 
Ef = k F /2m, where /cf = (37r 2 no) 1//3 is a convenient inverse length scale corresponding to a density 
no = 3 x 10 13 cm~ 3 as it is typical for the MIT 6 Li experiment [1]. (d) The effective scattering lengths 
as a function of the field detuning. The solid line is for the bare molecules (see the definition in the 
text) while the dashed one is for the atoms. 



checked that the dependence of the effective scattering length on the field detuning is in 
agreement with a different calculation in Ref. [16] under the two-channel model (we can only 
apply the two-channel model in this limit because of a large closed channel population [4, 17]). 
Experimentally, the scattering length between the dressed molecules can be measured from the 
collective excitations of the trapped Fermi gas [18] or from the in-trap radius of the condensate 
cloud [19]; while it is difficult to measure the scattering length between the bare molecules. 

The simplified equation (9) determines the distribution of the bosonic molecules. This 
solution, combined with Eq. (7), then fixes the distribution of the fermionic atoms. As a 
simple illustrating example, we use them to solve the fermi condensate shape function in 
a harmonic trap on the BEC side of the resonance. We take the values of the parameters 
corresponding to 6 Li, and assume a total of N ~ 10 5 atoms trapped in a time-independent 
potential V^r) = |mwr 2 with uj/2tt ~ 100Hz, as is typical in the experiments [1]. Figures 2(a) 
and 2(b) show the condensate shape functions in two different regions with the magnetic field 
detunings B — Bo given respectively by — 268G and — 107G. The first detuning corresponds 
to a point deep in the BEC region with (fc^a s ) _1 ~ 11, where a s is the atom scattering length 
at that detuning and hp 1 is a convenient length unit defined in the caption of Fig. 1. The 
second one corresponds to the onset of the bosonic region with the atom chemical potential 
[i ~ and (fcpds) -1 ~ 0.8. We have shown in Fig. 2 the number distributions for the bare 
molecules and the fermi condensate. One can see that these distribution functions are smooth 
in space, without the artificial cutoff at the edge of the trap as in the LDA result. The closed 
channel population (the total bare molecule fraction) is calculated to be about 33% and 0.1% 
respectively for these two detunings. 

An important goal of calculation of this simple example is to self consistently check the 
approximations made in our derivation. First, to derive the basic equations (5,6), we have 
assumed that the order parameter should be slowly varying over the size of the Cooper pairs. 
From Fig. (2), we see that the characteristic length for the variation of the order parameter 
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Fig. 2 - The radial distributions of of the bare-molecule density (nt), the fermi-atom density (n/), 
and the total density (ntot) in a harmonic trap with the field detuning B — Bo given by — 268G (a), 
and — 107G (b and c), respectively. All the densities are in the unit of no (see the value in the caption 
of Fig. 1) 



is typically of lOOfc^ 1 , while the size of the Cooper pair is well below kp 1 at these detun- 
ings [17]. Therefore, this approximation should be well satisfied for typical experiments. 
Second, from the basic equations (5,6) to the simplified equation (9), we have used the ap- 
proximation V e ff(r,t) ~ VhJr, t). To check the validity of this approximation, we calculate 
the effective potentials V e ff (r) and (r) (time- independent in this case) with our solutions 
of 0f, (r) and /k(r) from the stationary harmonic trap, and the results are shown in Fig. 3. 



It is clear that the difference 



K//(r)-T///(r) 



is small compared with the magnitude of 



|Ve//(r)| for different values of k when the atom chemical potential [i < 0, which justifies 
the approximation Vhf(r, t) ~ V e ff(r,t) in that region. One can also see that the relative 

error V e //(r) — V*,f(r) /\V e ff(r)\ goes up significantly (from roughly 1CP 4 to 1CP 1 ) when 

one goes from the field detuning — 268G to — 107G. If one goes further to the resonance point, 
this approximation eventually breaks down, and one needs to use the basic equations (5,6) 
instead of the reduced equation (9). 

In summary, we have derived a set of dynamical mean-field equations for evolution of 
strongly interacting fermionic atoms in any slowly varying potential traps, and discussed 
methods to solve these equations. We show that on the BEC side of the Feshbach resonance, 
this set of equations are reduced to a generalized dynamical GP equation. As an illustrating 
example, we solve the reduced equations in the case of a harmonic trap, which self-consistently 
verifies the approximations made in our derivation. 

This work was supported by the NSF award (0431476), the ARDA under ARO contracts, 
and the A. P. Sloan Fellowship. 



REFERENCES 



[1] C.A. Regal, M. Greiner and D.S. Jin, Phys. Rev. Lett. 92, 040403 (2004); M.W. Zwierlein et 
ai, Phys. Rev. Lett. 92, 120403 (2004); C. Chin et al, Science 305, 1128 (2004); J. Kinast et 



W. YI and L.-M. DUAN: DYNAMICAL MEAN-FIELD EQUATIONS FOR TRAPPED FERMIONS 7 



20 
10 




9 
6 
3 

-3 



x 10 



x 10 



(a) 






- V effi 




50 


100 150 




V 






-4 




x 10 




(c) 


- V e ff 









(b) S v k 

eff 



x 10 



(d) 8V k 



100 200 100 200 

V V 



-V e ff and the potential difference SV^ff 



V P 



eff 



eff 



Fig. 3 - The shapes of the effective potential 
(both in the unit of Ef) at the field detuning B — Bo — —26867 (a,b) and — 10767 (c,d), respectively. 
Figs, (a) and (d) each have twenty curves corresponding to the wave vector |k| varies from to lOfcf. 



al, Science 307, 1296 (2005); M.W. Zwierlein et al, Nature 435, 1047 (2005). G.B. Partridge, 

et al, Phys. Rev. Lett. 95, 020404 (2005). 
[2] D.M. Eagles, Phys. Rev. 186, 456 (1969); A.J. Leggett, in Modern Trends in the Theory of 

Condensed Matter ( Springer- Ver lag, Berlin, 1980), pp. 13-27; P. Nozieres and S. Schmitt-Rink, 

J. Low Temp. Phys. 59, 195 (1985). 
[3] M. Holland et al, Phys. Rev. Lett. 87, 120406 (2001); Y. Ohashi and A. Griffin, Phys. Rev. 

Lett. 89, 130402 (2002); 
[4] For a review, see Q.J. Chen, et al, Phys. Rep. 412, 1 (2005). 
[5] J. Kinnunen, M. Rodriguez, P. Torma, Science 305, 1131 (2004). 
[6] C. Caroli, P. de Gennes and J. Matricon, P hys. Lett. 9, 307 (1964). 

[7] Y. Ohashi, A. Gri ffin, |cond-mat/04102~2 preprint, (2004); C.-C. Chien, Y. H., Q.J. Chen, 
K. Levin, cond-mat/0510647 preprint, (2005); R. Sensarma, M. Randeria, T. L. Ho, 
|cond-mat / 05 1 076 1 1 preprint , (2005 ) . 
[8] J.R. Schrieffer, in Theory of Superconductivity (W.A. Benjamin, New York, 1964). 
[9] P. Pieri and G.C. Strinati, Phys. Rev. Lett. 91, 030401 (2003); A. Perali, P. Pieri, and G.C. 
Strinati, Phys. Rev. A 68, 031601 (R) (2003). 
[10] For a review, see Y. Castin, cond-mat/0105058 preprint, (2001) 
[11] M. Urban and P. Schuck, Phys. Rev. A 73, 013621 (2006). 
[12] E. Timmermans, et al, Phys. Rep. 315, 199 (1999). 
[13] L.-M. Duan, Phys. Rev. Lett. 95, 243202 (2005). 
[14] W. Yi, L.-M. Duan, Phys. Rev. A 73, 013609 (2006). 

[15] D.S. Petrov, C. Salomon, and G.V. Shlyapnikov, Phys. Rev. Lett. 93, 090404 (2004); P. Pieri, 

G.C. Strinati, |cond-mat /030742T] preprint , (2003). 

[16] J. Stajic, Q.J. Chen, and K. Levin, Phys. Rev. A 71, 033601 (2005). 

[17] R.B. Diener and T.-L. Ho, cond-mat/0405174 preprint, (2004); cond-mat/0404517 preprint, 

(2004). 

[18] M. Bartenstein, et al, Phys. Rev. Lett. 92, 203201 (2004). 
[19] T. Bourdel, et al, Phys. Rev. Lett. 93, 050401 (2004). 



